Defining the input and output files
library(knitr)
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.3.2
# Function to determine the number of rows until fixation is reached
time_to_fixation <- function(data, threshold = 0.9) {
# Find the index of the first row where the allele frequency is 0.9 or higher
fixation_index <- which(data$allele_frequency >= threshold)[1]
# Return the number of rows until fixation is reached
return(fixation_index - 1)
}
| Selection_coefficient | Fixation_probability | Avg_Fixation_time | Min_fixation_time | Max_fixation_time |
|---|---|---|---|---|
| s=0.1 | 0.4 | 35.55 | 29 | 39 |
| s=0.2 | 9.0 | 33.95 | 24 | 39 |
| s=0.3 | 27.8 | 26.65 | 19 | 34 |
| s=0.4 | 40.0 | 22.30 | 13 | 38 |
| s=0.5 | 51.3 | 14.45 | 10 | 23 |
| s=0.6 | 37.0 | 12.55 | 10 | 18 |
| s=0.7 | 46.5 | 9.95 | 6 | 16 |
| s=0.8 | 42.6 | 7.60 | 5 | 11 |
| Selection_coefficient | Avg_Length_Mb | Avg_window_freq | Min_freq | Max_freq | Avg_freq_variant_100k_window | |
|---|---|---|---|---|---|---|
| 1 | s=0.1 | 4.425 | 56.0 | 12.9 | 100 | 57.1 |
| 4 | s=0.4 | 4.575 | 75.9 | 18.6 | 100 | 76.5 |
| 2 | s=0.2 | 4.770 | 57.8 | 7.1 | 100 | 58.7 |
| 3 | s=0.3 | 4.785 | 71.7 | 20.0 | 100 | 72.7 |
| 8 | s=0.8 | 5.120 | 81.3 | 15.7 | 100 | 81.5 |
| 7 | s=0.7 | 5.210 | 84.8 | 31.4 | 100 | 85.6 |
| 5 | s=0.5 | 5.385 | 83.9 | 22.9 | 100 | 85.1 |
| 6 | s=0.6 | 5.485 | 78.0 | 10.0 | 100 | 78.5 |
Confidence level <=> konfidensgrad
# Function to calculate standard error and its confidence interval
standard_error_confidence_interval_fun <- function(observed_data, confidence_level = 0.95) {
# Calculate standard error
n <- length(observed_data)
standard_deviation <- sd(observed_data)
standard_error <- standard_deviation / sqrt(n - 1)
# Calculate confidence interval based on standard error
alpha <- (1 - confidence_level) / 2
margin_of_error <- qnorm(1 - alpha) * standard_error
mean_estimate <- mean(observed_data)
# Calculate the percentiles for the confidence interval
confidence_interval_lower_bound <- mean_estimate - margin_of_error # 2.5th percentile (2σ)
confidence_interval_upper_bound <- mean_estimate + margin_of_error # 97.5th percentile (2σ)
# Return confidence interval
return(c(confidence_interval_lower_bound, confidence_interval_upper_bound))
}
# # Function to calculate bootstrap confidence intervals
# standard_error_confidence_interval_fun <- function(observed_data, n_bootstraps = 1000, confidence_level = 0.95) {
#
# # Function to calculate the mean for each bootstrap sample
# compute_bootstrap_mean_fun <- function(observed_data_urn) {
# bootstrap_dataset <- sample(observed_data_urn, replace = TRUE)
# bootstrap_estimate <- mean(bootstrap_dataset)
# return(bootstrap_estimate)
# }
#
# # Perform bootstrap resampling
# bootstrap_point_estimates <- replicate(n_bootstraps, compute_bootstrap_mean_fun(observed_data))
#
# # Calculate the percentiles for the confidence interval
# alpha <- (1 - confidence_level) / 2
# confidence_interval_lower_bound <- quantile(bootstrap_point_estimates, alpha) # 2.5th percentile
# confidence_interval_upper_bound <- quantile(bootstrap_point_estimates, 1 - alpha) # 97.5th percentile
#
# # Return the confidence interval
# return(c(confidence_interval_lower_bound, confidence_interval_upper_bound))
# }
## ROH-hotspot selection testing results://n
| Model | Avg_ROH_hotspot_threshold |
|---|---|
| Neutral | 30.5 |
| s=0.1 | 32.2 |
| s=0.2 | 35.1 |
| s=0.3 | 37.3 |
| s=0.4 | 37.9 |
| s=0.6 | 41.4 |
| s=0.7 | 41.8 |
| s=0.5 | 42.6 |
| s=0.8 | 47.6 |
| Empirical | 75.0 |
## Overall Average Avg_F_ROH for german_shepherd : 0.2753012 //n
## Point estimate F_ROH across all 20 simulations: 0.25007 //n
## [1] "Bootstrap 95% Confidence Interval: [0.23473890934696, 0.265401133510183]"
## Point estimate F_ROH across all 20 simulations for selection_model_s01_chr3 : 0.2734885 //n[1] "Bootstrap 95% Confidence Interval: [0.259682877404132, 0.287294222595868]"
## Point estimate F_ROH across all 20 simulations for selection_model_s02_chr3 : 0.290046 //n[1] "Bootstrap 95% Confidence Interval: [0.268404249001561, 0.311687808141296]"
## Point estimate F_ROH across all 20 simulations for selection_model_s03_chr3 : 0.3161388 //n[1] "Bootstrap 95% Confidence Interval: [0.294101862625108, 0.338175708803464]"
## Point estimate F_ROH across all 20 simulations for selection_model_s04_chr3 : 0.3145064 //n[1] "Bootstrap 95% Confidence Interval: [0.286793720981368, 0.342219150447203]"
## Point estimate F_ROH across all 20 simulations for selection_model_s05_chr3 : 0.3394231 //n[1] "Bootstrap 95% Confidence Interval: [0.313208544370823, 0.36563759848632]"
## Point estimate F_ROH across all 20 simulations for selection_model_s06_chr3 : 0.348661 //n[1] "Bootstrap 95% Confidence Interval: [0.319210697417095, 0.378111345440048]"
## Point estimate F_ROH across all 20 simulations for selection_model_s07_chr3 : 0.3598673 //n[1] "Bootstrap 95% Confidence Interval: [0.337347371590621, 0.382387271266521]"
## Point estimate F_ROH across all 20 simulations for selection_model_s08_chr3 : 0.3911701 //n[1] "Bootstrap 95% Confidence Interval: [0.350246216491064, 0.432094012080364]"
| Model | F_ROH | Lower_CI | Upper_CI |
|---|---|---|---|
| Neutral | 0.250 | 0.235 | 0.265 |
| s=0.1 | 0.273 | 0.260 | 0.287 |
| Empirical | 0.275 | NA | NA |
| s=0.2 | 0.290 | 0.268 | 0.312 |
| s=0.4 | 0.315 | 0.287 | 0.342 |
| s=0.3 | 0.316 | 0.294 | 0.338 |
| s=0.5 | 0.339 | 0.313 | 0.366 |
| s=0.6 | 0.349 | 0.319 | 0.378 |
| s=0.7 | 0.360 | 0.337 | 0.382 |
| s=0.8 | 0.391 | 0.350 | 0.432 |
## Models where empirical F_ROH is within CI:
## [1] "s=0.1" "s=0.2"
##
## Models where empirical F_ROH is outside CI:
## [1] "s=0.3" "s=0.4" "s=0.5" "s=0.6" "s=0.7" "s=0.8" "Neutral"
## Uncommented because change of analysis
## Point estimate H_e across all 20 simulations for s01_chr3 : 0.05929176 //n[1] "Bootstrap 95% Confidence Interval: [0.0433130104967301, 0.0752705108634813]"
## Point estimate H_e across all 20 simulations for s02_chr3 : 0.05801723 //n[1] "Bootstrap 95% Confidence Interval: [0.0419902747617942, 0.0740441924354399]"
## Point estimate H_e across all 20 simulations for s03_chr3 : 0.03775639 //n[1] "Bootstrap 95% Confidence Interval: [0.0220355491960428, 0.0534772391642509]"
## Point estimate H_e across all 20 simulations for s04_chr3 : 0.03080788 //n[1] "Bootstrap 95% Confidence Interval: [0.0148546775692621, 0.0467610731950587]"
## Point estimate H_e across all 20 simulations for s05_chr3 : 0.02331459 //n[1] "Bootstrap 95% Confidence Interval: [0.00885300142719239, 0.0377761805510335]"
## Point estimate H_e across all 20 simulations for s06_chr3 : 0.03246187 //n[1] "Bootstrap 95% Confidence Interval: [0.0140181073288633, 0.0509056361685207]"
## Point estimate H_e across all 20 simulations for s07_chr3 : 0.01976834 //n[1] "Bootstrap 95% Confidence Interval: [0.00543074078237546, 0.0341059426184369]"
## Point estimate H_e across all 20 simulations for s08_chr3 : 0.02504322 //n[1] "Bootstrap 95% Confidence Interval: [0.0104102212736186, 0.0396762237506187]"
| Model | H_e_5th_percentile | Lower_CI | Upper_CI |
|---|---|---|---|
| Neutral | 0.12797 | 0.12039 | 0.13555 |
| s08_chr3 | 0.12831 | 0.12031 | 0.13632 |
| s01_chr3 | 0.12840 | 0.12246 | 0.13433 |
| s02_chr3 | 0.13008 | 0.12140 | 0.13876 |
| s04_chr3 | 0.13092 | 0.12240 | 0.13943 |
| s03_chr3 | 0.13443 | 0.12843 | 0.14043 |
| s06_chr3 | 0.13551 | 0.12673 | 0.14429 |
| s07_chr3 | 0.13766 | 0.12616 | 0.14916 |
| s05_chr3 | 0.14148 | 0.13333 | 0.14964 |
| Empirical | NA | NA | NA |
##
## ROH-hotspot threshold comparison between the different datasets
| Model | Avg_ROH_hotspot_threshold |
|---|---|
| Neutral | 30.5 |
| s=0.1 | 32.2 |
| s=0.2 | 35.1 |
| s=0.3 | 37.3 |
| s=0.4 | 37.9 |
| s=0.6 | 41.4 |
| s=0.7 | 41.8 |
| s=0.5 | 42.6 |
| s=0.8 | 47.6 |
| Empirical | 75.0 |
| Model | F_ROH | Lower_CI | Upper_CI |
|---|---|---|---|
| Neutral | 0.250 | 0.235 | 0.265 |
| s=0.1 | 0.273 | 0.260 | 0.287 |
| Empirical | 0.275 | NA | NA |
| s=0.2 | 0.290 | 0.268 | 0.312 |
| s=0.4 | 0.315 | 0.287 | 0.342 |
| s=0.3 | 0.316 | 0.294 | 0.338 |
| s=0.5 | 0.339 | 0.313 | 0.366 |
| s=0.6 | 0.349 | 0.319 | 0.378 |
| s=0.7 | 0.360 | 0.337 | 0.382 |
| s=0.8 | 0.391 | 0.350 | 0.432 |
## [1] "Selection test results"
## [1] "ROH-hotspot windows with an mean H_e Value lower or equal to the lower confidence interval of the fifth percentile of the neutral model are classified as being under selection"
## [1] "5th percentile of the neutral model is: 0.1203882"
| Name | Under_selection | Window_based_Average_H_e |
|---|---|---|
| Hotspot_chr3_window_1 | No | 0.12799 |
| Hotspot_chr3_window_3 | No | 0.13260 |
| Hotspot_chr17_window_2 | No | 0.14866 |
| Hotspot_chr3_window_2 | No | 0.15005 |
| Hotspot_chr17_window_1 | No | 0.15685 |
| Hotspot_chr19_window_1 | No | 0.17061 |
## [1] "ROH-hotspots under selection:"
| Name | Under_selection | Window_based_Average_H_e |
|---|
| Model | H_e | Lower_CI | Upper_CI | Under_Selection |
|---|---|---|---|---|
| s=0.7 | 0.01977 | 0.00543 | 0.03411 | Yes |
| s=0.5 | 0.02331 | 0.00885 | 0.03778 | Yes |
| s=0.8 | 0.02504 | 0.01041 | 0.03968 | Yes |
| s=0.4 | 0.03081 | 0.01485 | 0.04676 | Yes |
| s=0.6 | 0.03246 | 0.01402 | 0.05091 | Yes |
| s=0.3 | 0.03776 | 0.02204 | 0.05348 | Yes |
| s=0.2 | 0.05802 | 0.04199 | 0.07404 | Yes |
| s=0.1 | 0.05929 | 0.04331 | 0.07527 | Yes |
| Neutral | 0.12797 | 0.12039 | 0.13555 | No |
*** Behöver bootstrap av 5th percentiles från neutral model
| Model | Lengths_Mb | Avg_ROH_Freq |
|---|---|---|
| Hotspot_chr3_window_1 | 10.900 | 81.3 |
| s=0.6 | 5.485 | 78.0 |
| s=0.5 | 5.385 | 83.9 |
| s=0.7 | 5.210 | 84.8 |
| s=0.8 | 5.120 | 81.3 |
| s=0.3 | 4.785 | 71.7 |
| s=0.2 | 4.770 | 57.8 |
| s=0.4 | 4.575 | 75.9 |
| s=0.1 | 4.425 | 56.0 |
| Hotspot_chr19_window_1 | 4.400 | 75.6 |
| Hotspot_chr3_window_2 | 3.200 | 76.3 |
| Hotspot_chr3_window_3 | 2.700 | 77.6 |
| Hotspot_chr17_window_1 | 2.000 | 76.4 |
| Hotspot_chr17_window_2 | 0.600 | 76.1 |
## FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## Models where empirical H_e is within CI for Hotspot_chr3_window_1 :
## character(0)
##
## Models where empirical H_e is outside CI for Hotspot_chr3_window_1 :
## [1] "s=0.1" "s=0.2" "s=0.3" "s=0.4" "s=0.5" "s=0.6" "s=0.7" "s=0.8"
## FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## Models where empirical H_e is within CI for Hotspot_chr3_window_3 :
## character(0)
##
## Models where empirical H_e is outside CI for Hotspot_chr3_window_3 :
## [1] "s=0.1" "s=0.2" "s=0.3" "s=0.4" "s=0.5" "s=0.6" "s=0.7" "s=0.8"
## FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## Models where empirical H_e is within CI for Hotspot_chr17_window_2 :
## character(0)
##
## Models where empirical H_e is outside CI for Hotspot_chr17_window_2 :
## [1] "s=0.1" "s=0.2" "s=0.3" "s=0.4" "s=0.5" "s=0.6" "s=0.7" "s=0.8"
## FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## Models where empirical H_e is within CI for Hotspot_chr3_window_2 :
## character(0)
##
## Models where empirical H_e is outside CI for Hotspot_chr3_window_2 :
## [1] "s=0.1" "s=0.2" "s=0.3" "s=0.4" "s=0.5" "s=0.6" "s=0.7" "s=0.8"
## FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## Models where empirical H_e is within CI for Hotspot_chr17_window_1 :
## character(0)
##
## Models where empirical H_e is outside CI for Hotspot_chr17_window_1 :
## [1] "s=0.1" "s=0.2" "s=0.3" "s=0.4" "s=0.5" "s=0.6" "s=0.7" "s=0.8"
## FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## Models where empirical H_e is within CI for Hotspot_chr19_window_1 :
## character(0)
##
## Models where empirical H_e is outside CI for Hotspot_chr19_window_1 :
## [1] "s=0.1" "s=0.2" "s=0.3" "s=0.4" "s=0.5" "s=0.6" "s=0.7" "s=0.8"